home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Dr. Windows 3
/
dr win3.zip
/
dr win3
/
PROGRAMR
/
GSRC208A.ZIP
/
GAUSS.C
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-24
|
11KB
|
352 lines
#include "copyleft.h"
/*
GEPASI - a simulator of metabolic pathways and other dynamical systems
Copyright (C) 1989, 1992 Pedro Mendes
*/
/*************************************/
/* */
/* reduction of stoicheiometry by */
/* the Gauss method with row and */
/* column switches */
/* */
/* MICROSOFT C 6.00 */
/* QuickC/WIN 1.0 */
/* ULTRIX cc */
/* GNU gcc */
/* */
/* (include here compilers that */
/* compiled GEPASI successfully) */
/* */
/*************************************/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "globals.h"
#include "globvar.h"
#define ALMOST_ZERO 1.0e-7
void rowsw( int r1, int r2, int c )
{
int j;
float dummy;
double dumb;
int dum;
for(j=0;j<c;j++)
{
dum = stoi[r1][j]; /* switch rows in stoi and */
stoi[r1][j] = stoi[r2][j];
stoi[r2][j] = dum;
dummy = rstoi[r1][j]; /* switch rows in rstoi & */
rstoi[r1][j] = rstoi[r2][j];
rstoi[r2][j] = dummy;
}
dumb = x[r1]; /* switch rows in x and */
x[r1] = x[r2];
x[r2] = dumb;
dumb = x0[r1]; /* switch rows in x0 and */
x0[r1] = x0[r2];
x0[r2] = dumb;
j = ur[r1]; /* switch rows in ur too */
ur[r1] = ur[r2];
ur[r2] = j;
/*strcpy( buff, metname[r1] );
strcpy( metname[r1], metname[r2] );
strcpy( metname[r2], buff );*/
}
void colsw( int c1, int c2, int r )
{
int j;
float dummy;
int dum;
for( j=0; j<r; j++ )
{
dum = stoi[j][c1]; /* switch cols in stoi and */
stoi[j][c1] = stoi[j][c2];
stoi[j][c2] = dum;
dummy = rstoi[j][c1]; /* switch cols in rstoi & */
rstoi[j][c1] = rstoi[j][c2];
rstoi[j][c2] = dummy;
}
j = uc[c1]; /* switch cols in uc too */
uc[c1] = uc[c2];
uc[c2] = j;
/* strcpy( buff, stepname[c1] );
strcpy( stepname[c1], stepname[c2] );
strcpy( stepname[c2], buff );*/
}
void rowsw_m( int r1, int r2, int c ) /* switch rows in ml */
{
int j;
float dummy;
for(j=0;j<c;j++)
{
dummy = ml[r1][j];
ml[r1][j] = ml[r2][j];
ml[r2][j] = dummy;
}
}
void colsw_m( int c1, int c2, int r ) /* switch cols in ml */
{
int j;
float dummy;
for(j=0;j<r;j++)
{
dummy = ml[j][c1];
ml[j][c1] = ml[j][c2];
ml[j][c2] = dummy;
}
}
int emptyrow( int rw, int nsteps)
{
int j, ct; /* ct counts entries != 0 */
for( ct=0, j=0; j<nsteps; j++)
if ( fabs( rstoi[rw][j] ) > ALMOST_ZERO ) ct++;
return ( ct );
}
void invml11( void )
{ /* as ml is lower triangul.*/
int i,j,k; /* inverting is simply to */
float acum; /* forward-substitute with */
/* the identity matrix */
for(i=0;i<nmetab; i++) ml[i][i] = 1 ;
for( j=0; j<indmet; j++)
for( k=0; k<indmet; k++)
{
for( acum=0, i=0; i<k; i++ )
acum += ml[k][i] * lm[i][j];
lm[k][j] = ( (k==j) ? (float) 1.0 : - acum ) / ml[k][k];
}
}
void calc_ld( void )
{
int i,j,k;
for( i=indmet; i<nmetab; i++ ) /* first calculate L0 */
for( j=0; j<indmet; j++ ) /* which is */
{
ld[i][j] = 0;
for( k=0; k<indmet; k++ ) /* -1 */
ld[i][j] += ml[i][k] * lm[k][j]; /* L0 = L21 * (L11) */
}
for( i=indmet; i<nmetab; i++ ) /* Now map the lin. dep. */
{
for( j=0; j<indmet; j++ ) /* which is */
ld[i][j] = -ld[i][j]; /* -L0 */
for( j=indmet; j<nmetab; j++ ) /* concateneted with the */
ld[i][j] = (i == j)
? (float) 1.0 : (float) 0.0; /* m0 -> m part of Id. */
}
}
void init_moiety( void )
{
int i,j;
for( i=indmet; i<nmetab; i++)
{
moiety[i] = (float) 0.0;
for( j=0; j<nmetab; j++ )
moiety[i] += ld[i][j] * x0[j];
}
}
void lindep( void )
{
invml11(); /* invert ml11 into lm */
calc_ld(); /* calculate L2 * -L1' */
}
int gauss( void )
{
int i, j, k, flag;
float m;
for( k=0; k<nsteps+1; k++)
{
for( flag=0, i=k; i<nmetab; i++ )
if( fabs(rstoi[i][k]) > ALMOST_ZERO )
{
flag = 1;
break; /* suitable divisor */
}
if ( flag )
{
if ( i != k )
{
rowsw( k, i, nsteps ); /* switch row k with i */
rowsw_m( k, i, nmetab );
}
}
else
{
do
{
for( j=k+1; j<nsteps; j++ )
if( fabs(rstoi[k][j]) > ALMOST_ZERO )
{
flag = 1;
break; /* suitable value */
}
if ( flag && ( j != k ) )
{
colsw( j, k, nmetab ); /* switch col j with k */
colsw_m( j, k, nmetab );
}
if( !flag )
{
for( i=0; i<nmetab; i++ ) /* get rid of empty rows */
{
if ( ! emptyrow( i, nsteps ) ) /* all entries in row = 0 */
{
for( j=i+1; j<nmetab; j++)
if ( emptyrow( j, nsteps ) ) /* j = first non-empty row */
{
rowsw( i, j, nsteps ); /* switch row i with row j */
rowsw_m( i, j, nmetab );
flag = 2;
break;
}
}
}
}
if ( !flag )
{
indmet = k; /* # of independent metabs */
return(-1); /* flag conservation & ret */
}
}
while( flag != 1 );
}
for( i=k+1; i<nmetab; i++)
{
if ( fabs(rstoi[k][k]) > ALMOST_ZERO )
{
m = rstoi[i][k] / rstoi[k][k]; /* calculate the divisor */
ml[i][k] = m; /* and keep its symmetric */
if( m != (float) 0.0 )
{
for( j=k; j<nsteps; j++ ) /* reduce another row */
{
rstoi[i][j] = rstoi[i][j] - m * rstoi[k][j];
if( fabs(rstoi[i][j]) <= ALMOST_ZERO ) /* then make it a true zero*/
rstoi[i][j] = (float) 0.0;
}
}
}
}
}
return 0; /* job done: no lin. dep.! */
}
void int_ord( void ) /* start by putting all internal */
{ /* metabolites as the first */
int i,j,k; /* m rows in stoi */
for( i=0, k=0; i<totmet; i++ )
if (intmet[i])
{
ur[k] = i; /* record the old number */
for( j=0; j<nsteps; j++ )
stoi[k][j] = stoiu[i][j];
k++;
}
nmetab = k; /* the no.of internal metabolites */
}
void ext_ord( void ) /* now put the declared exernal */
{ /* metabolites in the last rows */
int i,j,k; /* of stoi */
k = nmetab;
for( i=0; i<totmet; i++ )
if (!intmet[i])
{
ur[k] = i; /* record the index */
for( j=0; j<nsteps; j++ )
stoi[k][j] = stoiu[i][j];
k++;
}
nextmet = totmet - nmetab; /* the no.of external metabolites */
for(i=0; i<nmetab; i++) intmet[i] = 1; /* set intmet